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ABSTRACT 


The waveguide mode tropospheric propagation effect prediction program, M- 
Layer, originally written by Naval Command Control and Ocean Surveillance Center, 
Research, Development, Test and Engineering Division (NRaD), is revised for 
greater accuracy, speed and stability. The accuracy improvement is achieved first by 
converting the extended complex number representation into the representation by 
the complex exponent then by re-writing the group of Airy function computation 
subroutines. This accuracy improvement makes it possible to implement a self- 
checking procedure for determining the proper method to evaluate the height gain 
function. Finally, a new mode locating algorithm is introduced which improves the 
efficiency of mode search and eliminates the looping problem observed. The revision 
has been documented and the new program source code has been delivered to 
NRaD. It is also recommended that the mode search protocol, not just the mode 
locating algorithm introduced in this revision, be completely revised for better 
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I, INTRODUCTION 


M-Layer is a FORTRAN program for computing the propagation factor of an 
electromagnetic (EM) wave in a stratified atmosphere. It is desirable to extend the 
capability of this program to include a layer of random medium representing the air- 
ocean interface. To achieve this goal, there are many basic theoretical problems 
which have to be answered. First of all, the effect of the earth curvature in this 
program is taken care of through the classical earth-flattening approximation [Ref. 1], 
but the result [Ref. 2] does not agree with the more recent diffraction theory of Fock 
[Ref. 3] near the surface of the earth. Then there is the question about the better 
method to model the atmospheric refractive index profile, either piecewise linear or 
quadratic, to be resolved by a new earth-flattening approximation under development 
at NPS. The new approximation will also determine the functions to be used for the 
representation of the EM fields in each layer through uniform asymptotic theories. 
Within some proper region, these new functions are expected to reduce to the Airy 
functions utilized by M-Layer. The evolutionary nature of this effort prompted this 
review to improve the inner workings of the M-Layer program. In particular, the 
subroutines to search for the modes and those for evaluating the Airy functions will 
remain as an important part of a program investigating questions about EM wave 


propagation by solving the related boundary value problem. 


It can never be overemphasized that a boundary value problem which includes 
a layer of random medium or some range dependent inhomogeneity, set up according 
to the Maxwell equations, will include backscattering in its solution. This is in sharp 
contrast to those numerical procedures based on the parabolic approximation to the 
wave equation for which the backscattering is completely ignored. 

In what follows, the M-Layer program and the reasons for replacing the 
extended complex numbers with their complex exponent representations are 
discussed, together with some other problems encountered and resolved during this 


investigation. 


A. M-LAYER 

In M-Layer, the index of refraction of the atmosphere is assumed to be height 
dependent and is soaceiatea with a continuous piecewise linear profile. The 
classical earth-flattening approximation Is utilized to allow the use of the cylindrical 
coordinate system while retaining the effect of the curvature of the earth. This is 
done simply by substituting the index of refraction with the modified index of 
refraction, which also has a piecewise linear profile [Ref. 1]. 

The source of the EM radiation is assumed to be either a vertical electric 
dipole or a vertical "magnetic dipole", with the latter providing an approximation to 
the radiation of a horizontal electric dipole. The dipole is located along the positive 
z-axis Of the cylindrical coordinate system while the origin is sitting on the ground. 


The x-y plane is the ‘flattened’ earth surface. After carrying out the Hankel 


transform along the radial direction, the resulting spectrum of the Hertzian dipole 
field within each layer of a linear segment of the modified refractive index profile is 
reduced to a linear combination of the Airy functions. Specifically, the layers are 
numbered to increase with height, with the first layer being the one above the 
ground. The spectrum of the Hertzian dipole field is proportional to the product of 
the values, at the transmitter height and at the receiver height respectively, of the 
height-gain function. At a height within the i-th layer, the height-gain function is 


given by [Ref. 4]: 


filp.2)=Bi(p)IA (pk, (9) +4 (4)) (1) 


where p is the radial component of the propagation vector and Is also the spectral 
variable of the Hankel transform; hence it is the same throughout all layers. It is a 
complex variable whose imaginary part represents the radial attenuation rate of the 
spectral component of the Hertzian dipole field. Under the classical earth-flattening 
approximation, the spectrum of the Hertzian dipole field contains a discrete portion 
and a branch cut. The discrete spectrum gives rise to the creeping wave modes 
diffracted by the earth surface and the dielectric waveguide modes supported by the 
layered atmosphere. The contribution from the branch cut is usually negligible, 
especially for the field in the shadow of the earth. The M-Layer program locates the 
discrete spectrum for modes having a radial attenuation rate below a predetermined 
value. Contributions from these modes determine the propagation factor of the 


Wave. 


The variable q; in the i-th layer is a dimensionless linear function of height z 
with the free space wavenumber k, the modified index of refraction m, at the lower 


boundary z = z;, the slope of the modified index of refraction a;/2 and p as 


3 k 2 2 
q;= a miva(e-s)-£5} (2) 
oF k? 


The height dependence of the field is given in terms of the functions k,(q;) and 


parameters: 


k,(q;), which are proportional to the Airy functions Ai(-q el? W/. 2) and Al(-q;) 
respectively. Of these two functions, at a height so large that q; is large and positive, 
k,q,;) represents a downward going wave and elf t/. k@) +k.(q;) represents an 
upward going wave. The coefficients. A ; andB,; are determined by the conditions on 
the continuity of the eseaan dipole field and its derivative across layer boundaries 
and by the normalization condition that the integral of the square of the height-gain 
function over all height equals unity. 

To fulfill the radiation condition, the highest layer is given the same refractive 
index as the free space above it and only the outgoing wave is allowed within this 
layer. Below the ‘flattened’ earth surface, the field is assumed to be a plane wave 
propagating downward. Hence, only the normalization factors are required in the 
highest layer and in the ground. By assigning B, to unity in the highest layer, all the 
coefficients A; and B; can be determined, according to the boundary conditions, to 


within a multiplicative factor for B;. This multiplicative factor is then deduced from 


the normalization condition. This procedure can also be carried out from the ground 
level up. That these coefficients can be computed either from the highest level down 
or from the lowest level up is a result of the fact that p belongs to the discrete 
spectrum of the Hertzian dipole field. Consequently, agreement between these two 
ways of evaluating the A, and B;, coefficients confirms that a mode has been located 


accurately. 


B. EXTENDED COMPLEX NUMBER REPRESENTATION 

The discrete spectrum of the Hertzian dipole field corresponds to the zeroes 
of the modal function which is a determinant whose elements consist of k ,(@q;) and 
k.(q;) at the layer boundaries. Numerically, the magnitude of this modal function 
causes Overflow and underflow problems as k,@;) or k»@;) becomes exponentially 
large or small for complex q; values. In the M-Layer program, to overcome this 
problem, a complex number is written as a scaled number, which is complex, 
multiplied by a scaling factor which is an integer power of e, the base of natural 
logarithm. This integer is chosen so that the greater of the absolute values of the 
real part and the imaginary part of the scaled number lies within e*!. A complex 
number written in this form is called an extended complex number. Multiplication 
of two extended complex numbers requires summing the two integer exponents in 
addition to carrying out the regular complex multiplication of the scaled numbers. 
Addition of two such numbers is achieved through the use of an addition subroutine: 


the larger scaling factor is factored out of both addends before they are combined. 


The scaling factor is adjusted after each addition and after a sequence of 
multiplications to make sure that the resulting scaled number is still within the 
desired range. Addition is troublesome when the two numbers to be added nearly 
cancel each other. Under this circumstance, the scaling factors of the two numbers 
are identical and both the real parts and the imaginary parts of the scaled numbers 
are almost equal with opposite signs. It is clear that the real part and the imaginary 
part of the sum lose their accuracies to different degrees; hence the phase angle may 
incur substantial error. To remedy this situation, interpolation procedures have to 
be devised. 

As two complex numbers come close to cancel each other, they must be out of 
phase by almost 180 degrees. By factoring out the square root of their product 
instead of the scale factor, the resulting addends become reciprocal to each other, 
both lying within an identical small angle to, and on the same side of, the imaginary 
axis. They are close to the unit circle, but one is on the inside and the other is on 
the outside. Taking out further a phase factor of 7/2 after writing the addends in 
their exponential forms, the exponents become small numbers for which a Taylor 
series expansion of the exponential function converges rapidly and can be used for 
interpolating the sum to achieve higher accuracy. Note that after the extra phase 
factor of 17/2 is removed from the addends, it is actually the difference of the 
resulting two reciprocals which is computed. This procedure effectively picks the 


direction on the complex plane along which the addends are almost opposing each 


other to carry out their cancellation. The resulting sum has a phase angle nearly 
perpendicular to this chosen direction. 

It is evident that the representation of a complex number by its complex 
exponent of base e provides better phase accuracy for addition. A one-to-one 
correspondence can be achieved by restricting the imaginary part of this exponent to 
within -7 and 7. This will be called the exponential representation or the complex 
exponent representation henceforth. It is convenient for multiplication: adding the 
complex exponents of the two factors will suffice. Conversion of the M-Layer 
program from the extended complex number to the complex exponent representation 


has been carried out. 


C. OTHER REVISIONS 

As better precision is achieved, problems with the mode search procedure and 
the evaluation of the A; and B; coefficients become severe. They are thoroughly 
investigated and resolved. For mode search, although the division of the region of 


interest into "contour rectangles" and further into square "meshes", and the search 


pattern to r the sides of a "contour rectangle" to find and follow "phase 
lines" int. ~ basic assumption of Shellman and Morfitt [Ref. 5] that 
both the r. zinary parts of the modal function are linear along every 
edge of a mi .< 1s completely abandoned. For the evaluation of the A; and 


B, coefficients, the "test for evanescence" conditions have been removed. A condition 


to determine whether to evaluate the coefficients from the ground level up or from 


the top level down has been fomulated and incorporated into the program. This 
accomplishment leads to the relaxation of mode locating accuracy requirements 
which, combined with the improved precision of the revised program, makes the first 
order Newton-Raphson iteration unnecessary. The specific changes in the program 
and the resulting gains in speed, accuracy and execution stability are discussed in the 
following chapters. Suggestions to completely revise the mode search protocol to do 
without the "contour rectangles" and to look for the modes according to their range 


attenuation rates are also provided. 


Il. PROGRAM REVISIONS 


M-Layer is structured into three parts: setup, mode search and propagation 
factor evaluation. The main input is the modified refractive index values at specified 
heights so that a piecewise linear profile can be constructed. If the mode locations 
for the particular profile are available from a previous run of the program, they can 
also be included in the input and the mode search procedures will be bypassed. The 
various ranges and transmitter and receiver heights for which propagation factors are 
desired are also specified. The subroutine WVGSTDIN is called to input the 
information from an ASCII data file. The program then computes the constants to 
be used for mode search and propagation factor evaluation. The mode search is 
performed with the subroutine FNDMOD. The MODSUM subroutine is then 
invoked to first compute the A; and B; coefficients as explained in the Introduction, 
then compute the propagation factor and the propagation loss. The complete 
program structure is given in Figures 1 and 2. There are several other subroutines 
which are not included in these and other figures, such as DHORIZ for computing 
the horizon distance between a transmitter and a receiver for reference purpose; 
CHKMOD, a maintenance routine for removing zeroes from reported mode 
locations by older versions of the program; or AO2H20, a routine to compute the 


atmospheric absorption coefficient due to oxygen and water vapor. They will not be 





Figure 1 Original M-layer subroutines structure. 





Figure 2 Original M-layer subroutines structure (continued). 
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discussed as they do not contribute directly to the main purpose of this program of 
locating the modes and computing the propagation factor. 

The program structure has been altered as shown in Figures 3 and 4. Since the 
A; and B; coefficients have to be evaluated only once, they are now obtained through 
a call to the subroutine ABCOEF directly from the main program right after the 
modes are located. Several subroutines are dropped in this revision for various 
reasons: The subroutines NORME and NORMRE are eliminated because they are 
no longer needed due to the change in complex number representation; the 
subroutines NOMSHX, FDFDTX and DXDETR are not used because the modes 
are now located with adequate precision without further iteration; the subroutine 
ADDxX is not listed separately because it is called only once and has been reduced 
to only a few lines which are placed where the subroutine is called in the original 
program. On the other hand, changes in the mode search algorithm require the 
addition of two new subroutines: SURFO is a modified and simpler version of SURF; 
ROOTS replaces QUAD. Due to the change in complex number representation, all 
subroutines listed below FNDMOD and MODSUM have been revised, including 
their input/output lists. But except for SURFO and ROOTS, the utilities of these 
subroutines are the same as those of the original ones. Descriptions of these 
subroutines can be found in the report by Yeoh [Ref. 4]. 

The most significant changes have been made in XCADD, XCDAIT and 


XCDAIG for adopting the complex exponent representation and improving 


Fe 





Figure 3 New M-lavyer subroutines structure. 


® © 
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Figure 4 New M-layer subroutines structure (continued). 
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computation speed and accuracy; in FZEROX, FINDFX, ROOTS and SURFO for 
stabilizing and simplifying the mode search algorithm; and in ABCOEF for 
implementing the criteria to determine the reliable manner for evaluating the A; and 
B; coefficients. These changes are discussed in the sections below. The source code 
listings of the completely new subroutines XCADD and ROOTS and the significantly 
revised subroutines FZEROX and ABCOEF, which are compiled with Microsoft 
FORTRAN version 5.00, are attached as Appendices A through D. Validation of 
the revised program has been carried out at 9.6 GHz for all the 21 profiles listed in 


Yeoh [Ref. 4]. 


A. ADDITION SUBROUTINE 

XCADD is the subroutine implementing the addition of complex numbers 
under the representation by their exponents. Given the double precision complex 
numbers z, and z, as the exponents of the addends, this subroutine returns the 
exponent of the sum. Since a double precision number has an accuracy of 53 bits, 
if the real parts of z, and z, differ by more than 53 bits, the exponent of their sum 
will simply be the one of the greater real part. When cancellation becomes Serious, 
the square root of the addends is factored out first. Then the four-term Taylor series 
expansions of the resulting reciprocals are summed. Since the leading term of the 
suin of the Taylor series is a good estimate of the sum of the reciprocals and the 


relative error of the four-term Taylor series sum is proportional to the fourth order 
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of this leading term, the threshold for invoking this interpolation procedure is set at 


the highest possible value of pet 


allowed under double precision. Experimenting 
with this procedure shows that this interpolation improves accuracy as long as the 


threshold is set at a number between 2774 and 2714. 


B. AIRY FUNCTION EVALUATION 

Similar to the original program, the evaluation of the Airy function adopted the 
algorithm prescribed by Schulten, et. al. [Ref. 6]. In the new program, changes are 
made to follow the advice of Schulten, et. al. concerning the region within which a 
Taylor series expansion, instead of the faster Gaussian quadrature, has to be used to 
achieve double precision accuracy. Other changes in implementing the algorithm are 


described below. 


1. XCDAIT 
Due to the similarity in their Taylor series coefficients, the Airy function 
and its derivative are evaluated within a single loop. The relative accuracy of the 


derivative of the Airy function is set at the double precision limit of ee 


2. XCDAIG 
Six term Gaussian quadrature ts used for evaluating the Airy function and 
its derivative outside the circle of radius 4.97 centered at (0.90, 2.80) on the complex 
plane. The use of four-term quadrature outside a radius of 15 from the origin 


suggested by Schulten, et. al. is not adopted. The six-term quadrature in this range 
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retains a higher accuracy while overall speed improvement by using both the four- 


term and the six-term quadrature appears to be minimal. 


C. MODE LOCATING 

As explained in the Introduction, the modes are located at the zeroes of the 
modal function. These zeroes are located on the upper complex q,; plane. Here qy, 
is the value of q, on the earth’s surface, which, according to Eq.(2) of Chapter I, is 
a linear function of p*. For a horizontally propagating mode, p/k is close to unity. 
The maximum range attenuation rate specified for the desired modes, which 
corresponds to a limit on the imaginary part of p, determines approximately the 
upper bound for the imaginary part of the q,, complex plane to be searched for 
modes. The Shellman and Morffit mode search procedure first divides the search 
region horizontally into "contour rectangles" each of which spans 160 meshes along 
the real q,, direction. A mesh is a square whose size is an adjustable parameter of 
the order 1074 at 9.6 GHz for most of the cases considered herein. This parameter 
is determined by the frequency and the slope of the modified index of reflection in 
the lowest layer of the profile. The search commences at the top left corner of the 
“contour rectangle” whose left edge has a real coordinate value close to the 
difference of the real parts of the q,, values, with the minimum modified index of 
refraction and the index near the surface substituted into Eq.(2) of Chapter I. After 


the search over the initial rectangle is completed, the program moves to search the 


sh) 


next rectangle until a specified maximum number of modes are found or a specified 
number of "contour rectangles" have been searched. 

The search for zeroes makes use of the fact that a real function changes sign 
when it crosses a simple zero. Since a zero of a complex valued function F(q) is 
where both its real part and imaginary part vanish, a necessary condition for a point 
q, to be a zero is that it is on the intersection of two curves defined by Im{F(q)} = 
0 and Re{F(q)} = 0. The program searches around a "contour rectangle" for a sign 
change in Im{F(q)} across an edge of a mesh bordering the side of the "contour 
rectangle” to determine that a line of Im{F(q)} = 0 has been encountered. The 
search then follows this line into the meshes within the "contour rectangle", checking 
each mesh to see if a curve Re{F(q)} = 0 enters the mesh under investigation. All 
these steps make use ine of the assumption that the zeroes of the modal function 
are simple. Once both the curve Im{F(q)} = 0 and the curve Re{F(q)} = 0 are 
determined to be present within a mesh, the location of their possible interception 
is estimated. An algorithm for this estimate is required. 

Shellman and Morffit [Ref. 5] introduced a further assumption that the 
functions Re{F(q)} and Im{F(q)} are both linear along the edges of a mesh. Based 
on this assumption, they try to estimate the locations where the curve Im{F(q)} = 
Q enters and leaves a mesh square, and the location of q, if a curve Re{F(q)} = 0 
also enters the same mesh. It is obvious that information about the locations where 


the curves enter and leave the mesh square is not essential. Furthermore, in the 18 
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m duct height case, the scheme causes the search path to loop around four 
contiguous meshes until the search is broken up by the limit on the number of 
meshes to be investigated. Replacing their technique requires major changes in the 
subroutines involved. A new subroutine ROOTS is provided to estimate the location 
of the intersection of the curves Im{F(q)} = 0 and Re{F(q)} = 0. These changes 
eliminate the looping problem. 

Another problem is encountered in the 40 m duct height case when a large 
number of zeroes are found in the lower half complex q,, plane. These zeroes 
appear to belong to the reflection coefficient on the wrong sheet of the branch cut 
and are not waveguide modes. This happens because the search region has been 
extended below the real q,, axis to avoid the singularity in SURF. The problem with 
this singularity should have been solved within SURF, especially because it occurs 
only when the derivative of the subroutine output variable gamma with respect to q,, 
is computed. Since this derivative is not needed during mode search, the extension 
of the search region to the negative q,, plane is unnecessary. A simplified routine, 
SURFO, is introduced which is exactly the same as SURF except that it does not 
evaluate the derivative of gamma. By using this subroutine instead of SURF, the 


search path in the revised program does not avoid the real and the imaginary axes. 


1. FNDMOD 
The search region is limited to the upper half q,, plane. All the modes 
found are ordered according to their range attenuation rates before those numbered 


beyond the maximum modes allowed are abandoned. 


i 


2. FZEROX 

Since the curve Im{F(q)} = 0 enters into a mesh square through an edge, 
the values of Im{F(q)} must change sign over the end points of either one or all 
three other edges. When there is only one other edge across which Im{F(q)} 
changes sign at its end points, it is the edge across which the curve Im{F(q)} = 0 
exits the mesh square. Ambiguity arises when all edges indicate a change of sign at 
their end points. When this occurs, a "right turn rule" is adopted which assumes that 
the curve exits the edge to the right of the one along which it enters the mesh 
square. Such a rule avoids the retracing of the search path when the mesh square 
is revisited as entering this same mesh square from the left side of an edge after 
exiting from its right side requires a crossing of the Im{F(q)} = 0 curve, which is 
prohibited under the simple zero assumption. On the other hand, the actual curve 
may have turned left and then returns to this mesh square, 1.e., following a "left turn 


rule." Under such a scenario, this wrong choice would have left a segment of the 
curve not searched. This difficulty has not been observed during testing. In fact, the 
ambiguous situation seldom occurs. Note also that, as remarked above, two lines of 


Im{F(q)} = 0 do not cross each other unless a higher order zero is present. Hence, 


only a "right turn rule" or a "left turn rule” for the curve to exit the mesh is allowed. 
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Exiting the opposite edge demands a pair of crossing Im{F(q)} = 0 curves within the 
mesh square. This violates the assumption that all zeroes are simple. Also note that, 
the possibility of vanishing Re{F(q)} or Im{F(q)} values at the corners of a mesh 


square is eliminated through a small adjustment in FINDFX. 


3. FINDFX 

Both the vertical shift away from the real q,, axis and the horizontal 
offset away from the imaginary axis are unnecessary and have been removed from 
this routine. Furthermore, as a result of converting to the complex exponent 
representation, the sine and cosine of the argument of the modal functions are 
examined for sign changes in FZEROX. This is implemented in FINDFX by 
including the cosine and sine values of the argument of the modal function in the 
output list. To avoid the indeterminate case when either the real or the imaginary 
part of the modal function becomes zero at any corner of a mesh square, the 
argument for computing the cosine and sine values is increased by 2~>> when this 
occurs. This is equivalent to a consistent small distortion of the particular corner of 
the mesh square. This will not cause any error in locating the zero because FINDFX 


still returns separately the unmodified exponent of the value of the modal function. 


4. ROOTS 
Assuming that the modal function is analytic within the mesh, this 


Subroutine utilizes the values of the modal function at the four corners of the mesh 
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square to determine the Taylor series expansion coefficients of the modal function 
to the third order. The roots of this cubic polynomial are then located using 
Cardan’s solution by radicals. If the higher order coefficients fall below machine 
resolution for a root within the mesh square, these coefficients are regarded as zero 
and the order of the polynomial is reduced and can be solved more expediently. If 
the function is determined to be constant over the mesh square, the center of the 


Square is taken as the root location. 


D. EVALUATING A; AND B; 

As discussed in the Introduction, the A; and B; coefficients can be evaluated 
either from the top level down or from the lowest level up. These two procedures 
are simply called "integration down" and "integration up", respectively, in the original 
documentation [Ref. 4]. The location of a mode has been called an eigenvalue. 
That the results of "Integration down" and "integration up" agree is a manifestation 
that the eigenvalue is located accurately. 

The subroutine ABCOEF evaluates the coefficients A; and B; for each mode. 
If the range attenuation rate for a mode is greater than 0.1 dB/km, the coefficients 
are evaluated from the lowest layer up. Otherwise, it is evaluated from the top layer 
down. It is obvious that such a rule must be implemented because the results of 
"Integration up" and "integration down" do not agree for many modes. Efforts are 


made to determine the cause of this discrepancy and to devise a means to resolve it. 
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Investigation reveals that inadequate precision in the location of the modes is 
one source of the problem. Since the B; coefficients depend on the A; coefficients, 
while the A; coefficients are obtained directly, only the A; coefficients need to be 
examined. The A; coefficients of the six modes of lowest range attenuation rates for 
all 21 profiles except the one without evaporation duct are computed using 
eigenvalues of different accuracy controlled by the first order Newton-Raphson 
iteration method. Table 1 shows the A; coefficient computed with the new program. 
They are arranged from the top layer down. In the i-th layer, the A; coefficient 
computed by "integration downward" depends only on A;, ; in the layer above while 
that computed by “integration upward" depends only on A, _, in the layer below. 
Hence in each layer, the coefficient obtained by "Integration downward" Is listed 
above that obtained by "integration upward". There are five sets of A; values listed, 
with the magnitudes given in powers of 10, and the phase given as a multiple of 7. 
They are obtained from eigenvalues of decreasing accuracy -the one used to compute 
the left most column being the most accurate. The first set is computed using an 
eigenvalue having a relative accuracy of 2~% | The second set uses an eigenvalue with 
a relative accuracy of 2~*°. The relative accuracy of the eigenvalue for the third set 
is 2~*©. For the fourth set, the first order Newton-Raphson iteration of the mode 
location is set at an absolute accuracy of 0.03 of the mesh size, same as that specified 
in the original program. The eigenvalue for the right most set is the mode location 


estimated by ROOTS without modification by the Newton-Raphson iteration. It is 


ZA 


TABLE 1. IMPROVING A, ACCURACY WITH EIGENVALUE 18 M DUCT 


mode 4 q-etgenvalue: 
eigenvalue difference: 
layer # Ai/down 
layer # Ai/up 
18 .0261 .6719 
18 -0261 .6719 
17 -.0625 .6368 
17 -.0625 .6368 
16 -0139 .7440 
16 .0139 .7440 
15 -1216 .6353 
15 -1216 .6353 
14 -0166 .5471 
14 -0166 .5471 
13 -.1565 .5310 
13 -.1565 .5310 
i2 -.3842 .5659 
12 -.3842 .5659 
11 -2.2002 -.8081 
11 -2.2002 -.8081 
10 -5.4648 .2423 
10 “5.4648 .2423 
9 -3.6974 -.6979 
9 -3.6978 -.6978 
8 -3459 -.7982 
8 -3459 -.7983 
7 -4098 .8794 
7 -4097 =.8793 
6 -3480 = .8161 
6 -3479 = .8160 
5 -2923 .8304 
5 -2922 .8303 
4 -2399 #.8619 
4 -2398 .8618 
3 -1831 .8910 
3 -1831 .8908 
2 -1300 8 .9149 
2 -1300 .9146 
1 -0586 .9335 
1 -0586 .9331 


- 18885 743251768030+00 
.000+00 .00D+00 -.49D-13 -.660-15 
Ai/down Ai/down 
Ai/up Ai/up 

.0261 -6719 -0261 .6719 
.0261 .6719 .0261 -6719 
-.0625 .6368 -.0625 .6368 
-.0625 .6368 -.0625 .6368 
.0139 .7440 -0139 .7440 
.0139 .7440 .0139 .7440 
-1216 3.6353 .1216 .6353 
1216 .6353 -1216 .6353 
.0166 .5471 -0166 .5471 
.0166 .5471 -0166 = .5471 
-.1565 .5310 -.1565 .5310 
-.1565 .5310 -.1565 .5310 
-.3842 .5659 -.3842 .5659 
-.3842 .5659 -.3842 .5659 
-2.2002 -.8081 -2.2002 -.8081 
-2.2002 -.8081 -2.2002 -.8081 
-5.4648 .2423 -5.4648 .2423 
-5.4648 .2423 -5.4648 .2423 
-3.6974 -.6979 -3.6974 -.6980 
-3.6978 -.6978 -3.6978 -.6978 
3459 -.7982 -3460 -.7982 
3459 -.7983 -3459 -.7983 
-4098 .8794 .4098 .8794 
-4097 ~=.8793 -4097 = .8793 
-3480 .8161 -3480 =. 8161 
-3479 = .8160 -3479 .8160 
-2923 .8304 .2923 ~=—.. 8304 
-2922 ~3.8303 -2922 +#.8303 
-2359 =. 8619 .2360 .8619 
-2358 .8618 .2358 .8618 

- 1831 .8910 .1832 ~=.8910 
.1831 = .8908 -1831 .8908 
-1300 .9149 -1301 9149 
-1300 .9146 -1300 .9146 
-0586 = .9335 .0588 .9335 
.0586 = .9331 -0586 .9331 
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- 10806787448105980-01 


-12D-11 -.12D-10 


Ai/down 
Ai/up 

.0261 += .6719 
0261 6719 
-.0625 .6368 
-.0625 .6368 
-0139 .7440 
-0139 .7440 
21216 ~26353 
-1216 .6353 
-0166 .5471 
-0166 .5471 
-.1565 .5310 
-.1565 .5310 
-.3842 .5659 
-.3842 .5659 
-2.2002 -.8081 
-2.2002 -.8081 
-5.4654 .2423 
*5.4648 .2423 
-3.6783 -.7012 
-3.6978 -.6978 
3482 -.7926 
-3459 -.7983 
-4136 .8836 
-4097 # .8793 
-3526 = .8205 
3479 ~=.8160 
2972 =. 83558 
2922 #8=.8303 
-2408 .8690 
2358 .8618 
-1878 .9003 

- 1831 -8908 
1342 8.9275 
.1300 .9146 
-0618 .9545 
-0586 = .9331 


-15D-06 .60D-0 
Ai/down 
Ai/up 
-0261 6719 
.0261 6719 
-.0625 .6368 
-.0625 .6368 
.0139 .7440 
-0139 =.7440 
-1216 .635 
-1216 .635 
-0166 .5471 
-0166 = .5471 
-.1565 =.5310 
-.1565 .5310 
-.3843 .5659 
-.3842 .5659 
-2.1909 -.8068 
-2.2002 -.8081 
-4.1810 -.2161 
-5.4647 .242 
-6.4611 -.2121 
-3.6977 -.6978 
-1.9078 -.9148 
-3459 -.798 
-1.0899 .5364 
-4097 ~=.879 
-.5879 .4005 
-3479  .8160 
-.3490 .3749 
.2922 ~=—. 830 
-.2058 .3731 
.2358 .8618 
-.1250 = .375 
- 1831 -8908 
-.0734 .3750 
-1300 = .9146 
-.0318 .3670 
-0586 = .9331 


clear that, for this mode, the difference between these two methods of computing the 
coefficients becomes negligible as the accuracy in mode location increases. For 
example, in the 8-th layer, the magnitude of A; computed by integrating downward 
changes from —1.9078 to 0.3482, to 0.3460 to 0.3459, which agrees with the result 
computed by integrating upward. The phase follows the same trend to an agreement 
within 0.0017. Table 2 shows a similar set of output, but the coefficients fail to agree 
even when the relative accuracy is increased to 2 40. Note that the actual difference 
in both the real part and the imaginary part of the two most accurate eigenvalues Is 
about 2~*8. Double precision accuracy appears to be insufficient for the coefficients 
computed with these two methods to agree for all modes. Some interesting features 
can be observed in both tables, which are present in all 120 sets of values computed. 
When disagreement is present in one set of A; coefficients, such as those in either 
Table 1 or Table 2, the change toward smaller differences with improving eigenvalue 
accuracy occurs mainly in one way of computation, but not both. For example, in 
Table 1, the values of “integration downward" improve with better eigenvalue 
accuracy, while those computed by "integrating upward" change little. In Table 2, the 
results of "integration downward" are the ones that are holding steady as the accuracy 
in eigenvalue improves. Furthermore, when disagreement occurs, the layer in which 
the A; coefficient has the smallest magnitude, i.e., the one having the most negative 
power of 10, divides the table into two parts. The results of two different ways of 
computation agree in the layers above this one if they disagree in those below it, and 


vise versa. No explanation will be attempted. Instead, practical rules are drawn up 
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TABLE 2. IMPROVING A, ACCURACY WITH EIGENVALUE 36 M DUCT 


eigenvalue difference: 


mode 3 
layer # 
layer # 
27 - 
27 

26 

26 - 
25 

25 -1 
24 

24 = 76 
23 -2 
23-14, 
22. -\2e 
22 4-23 
Alle S> 
21 = -44 
20 -131 
20 -129 
19 -25 
19 -25 
18 -13 
18 -13 
17 -7 
17 -7 
16 -3 
16 ah 
15 -2 
15 -2 
14 -1 
14 -1 
13 -1 
13 -1 
12 : 
12 2 
11 - 
11 - 
10 - 
10 - 


q-eigenvalue: 


Ai /down 
Ai/up 

.0009 .6663 
-2595" «wroee 
.0007 .6678 
0111 - 3659 
0022 .6657 
8851 3913 
.0001 . 6809 
-4914  .6081 
-9495 = .5951 
5340 .7973 
1956 .9278 
.5827 -.9407 
-2395 -.2502 
-4517 -.8199 
.3304 -.9570 
-0146 -.2961 
-6088 -.9230 
-6090 -.9228 
-6970 .6510 
-6970 .6510 
-0384 .4145 
-0384 .4145 
3146 3.2991 
3146 3.2991 
-3132 =. 2632 
3132 -2632 
-5669 .2415 
-5669 = .2415 
-0838 .2352 
-0838 .2352 
6983 .2432 
6983 .2432 
3794 =. 2712 
tot oGemeeerie 
-0102 .3619 
-0102 .3619 


- 14796229405 720070 - 02 


-3148000164781392D+00 
-380-14 .360-14 -380-14 .360-14 
Ai/down Ai /down 
Ai/up Ai/up 
-.0009 .6663 -.0009 .6663 
.2353 =. 7582 2353 =. 7582 
.0007 .6678 .0007 .6678 
-.0111 .3659 -.0111 .3659 
.0022 .6657 .0022 .6657 
-1.8851 .3913 -1.8851 .3913 
.0001 6809 .0001 .6809 
-7.4914 .6081 -7.4914  .6081 
-2.9495  .5951 -2.9495  .5951 
-14.5340 .7973 -14.5340 .7973 
-12.1956 .9278 -12.1956 .9278 
-23.5827 -.9406 -23.5827 -.9406 
-35.2395 -.2502 -35.2395 -.2502 
-45.8691 .8599 -45.8691 .8599 
-131.3304 -.9570 -131.3304 -.9570 
-127.6070 .0248 -127.6070 .0248 
-25.6088 -.9230 -25.6088 -.9230 
-25.6184 -.9241 -25.6184 -.9241 
-13.6970 .6510 -13.6970 .6510 
-13.6970 .6510 -13.6970 .6510 
-7.0384 .4145 -7.0384 .4145 
-7.0384 .4145 -7.0384 .4145 
-3.3146 .2991 -3.3146 3.2991 
-3.3146 = .2991 -3.3146 .2991 
-2.3132 .2632 -2.3132 .2632 
-2.3132 .2632 -2.3132 .2632 
-1.5669 .2415 -1.5669 .2415 
-1.5669 .2415 -1.5669 .2415 
-1.0838 .2352 -1.0838 .2352 
-1.0838 .2352 -1.0838 .2352 
-.6983 .2432 -.6983 .2432 
-.6983 .2432 -.6983 .2432 
-.3754 ~=.2712 -.3754 = .2712 
-.3754 .2712 -.3754 = .2712 
-.0102 .3619 -.0102 .3619 
-.0102 .3619 -.0102 .3619 
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-.160-09 -.22D-09 


Ai/down 
Ai/up 
.0009 .6663 
2353 3=. 7582 
-0007 £.6678 
0111 - 3659 
0022 £.6657 
8851 3913 
.0001 6809 
-4914 .6081 
9495 = .5951 
5340 .7973 
.1956 .9278 
-5827 -.9406 
.2395 -.2502 
-4590 -.1252 
.3304 -.9570 
9124 -.9081 
-6088 -.9230 
5644 -.8054 
-6970 .6510 
-0618 .7675 
-0384 = .4145 
0308 .4179 
-3146 = =.2991 
-3146 3=.2991 
3132 8 38=.2632 
-3132 =.2632 
-5669 = .2415 
5669 = .2415 
-0838 .2352 
.0838 .2352 
-6983 .2432 
6983 .2432 
3794 -2fle 
3794 -2712 
.0102 .3619 
.0102 .3619 


t 


-53D-07 .280-0 

Ai/down 
Ai/up 

-.0009 .666 
2353 =. 7582 
.0007 .6678 
-.0111 3659 
-0022 #£.665 

-1.8852 .391 
.0001 .6809 
“7.4914 -6081 
-2.9495  .5951 

-14.5340 .79 
“12.1956 .9278 
-23.5827 -.9406 
-35.2396 -.2501 
-47.4594 -.1251 
131.3307 -.9569 
120.9305 .8279 
-25.6088 -.9230 
-20.2166 .7391 
-13.6970 .6510 
-10.8148 -3440 
-7.0384 .4145 
-6.3129  .1800 
-3.3146 3.2991 
-3.3116 .2970 
-2.3132 -2632) 
-2.3127 .2629 
-1.5669 .2415 
-1.5668 .2415 
-1.0838 .2352 
-1.0838 .2352 
-.6983 .2432 
-.6983 .2432 
-.3754 .2712 
-.3754 <27ie 
-.0102 3619 
-.0102 - 3619 


to take advantage of these facts. In Table 1, the process of "integration upward" goes 
through the troublesome 10-th layer and produces results which agree with the results 
of downward integration before the downward process goes through the 10-th layer. 
On the other hand, the downward integration is tripped up going across the 10-th 
layer and produces results which fail to agree with the results from upward 
integration. It is clear that the results from upward integration are the correct ones. 
This conclusion is further supported by the fact that improving the accuracy of the 
eigenvalue does not change significantly the results of upward integration. Similar 
argument leads to the conclusion that in Table 2, the results of downward integration 
are the correct values. 

It can be concluded from the above observations that one of the methods of 
computing the A, coefficients converges to the correct value much faster than the 
other. It is also found that this method of faster convergence is always able to arrive 
at the correct values for A; for all the cases under investigation. 

Table 3 lists the statistics of the method of integration which yields the correct 
A, coefficients for each of the 120 modes investigated. The differences in magnitudes 
and phases in the lowest layer and in the layer below the highest are also listed. 
Since for most of the cases, when disagreement in A, values occurs, the correct 
integration is upward -this is used as the default. To decide that downward 
integration should be utilized, the following steps are taken: The first A; value of 
downward integration is computed and compared to the value from upward 


integration. If the magnitudes in dB disagree by less than 0.02 dB, their phases will 
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be checked. If the phases differ by less than 10”, the agreement is deemed 
acceptable and the A; and B; coefficients computed from the lowest layer up are 
used. Otherwise, the coefficients are re-evaluated again from the highest layer down. 

Once the correct method of evaluating the A; and B; coefficients is used, the 
accuracy of the mode location becomes less critical. For all the cases investigated, 
the A; coefficients obtained from mode locations estimated with or without the 
Newton-Raphson first order iteration differ only by 0.06 dB in magnitude and 
0.00137 in phase at most. In fact, few cases show differences more than 0.002 dB 
and 0.00017. The Newton-Raphson iteration is not needed. Hence the subroutines 


NOMSHX, FDFDTX and DXDETR are removed. 
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TABLE 3. STATISTICS oS STATISTICS FOR EVALUATING A. COEFFICIENT | 


A|A,| (dB) Aarg(A,)/x 
Duct Mode # Evaluating Method 


height 





Layer Layer 
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TABLE 3. CONTINUED 2. 
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13.456 0.0311 


1.014 0.2347 
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III. CONCLUSIONS AND RECOMMENDATIONS 


A. Performance 

This revision of M-Layer converts the extended complex number representation 
of an exponentially large or small number into the direct representation by its 
complex exponent. The accuracy of the computation has been improved in two ways: 
First, an interpolation algorithm has been devised when severe cancellation of the 
addends is detected. Secondly, accuracy for the evaluation of the Airy function has 
been improved, not just by summing the Taylor series to double precision resolution 
and by adopting six-term Gaussian quadrature, but also by expanding the region 
within which the more expedient Gaussian quadrature is excluded in favor of the 
more accurate, but time-consuming, Taylor series summation. The improvement in 
accuracy 1S most easily seen from Table 1. 

As discussed in the Introduction, evaluating the A; and B; coefficients either 
from the lowest layer up (integration up), or from the top layer down (integration 
down), must result in the same values. This property provides a consistency check 
for the accuracy of the computation. For the six modes of lowest range attenuation 
rates of the 20 profiles of different duct heights, Table 1 lists the maximum 
difference for each mode which shows a discrepancy between these two methods of 


evaluating the A, coefficients. For each profile, the maximum value in magnitude 
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TABLE 1. MAXIMUM DIFFERENCE IN A, COEFFICIENT BETWEEN 
_ INTEGRATION UP AND DOWN 


Difference in A coefficient 
Duct ag —_— in Phase ae over 
| height ea re 1 
new 


ioe Oe ns i 
2 oe es ee ee 
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TAB LE 1. CONTINUED 


Difference in A coefficient 


Magnitude difference in Phase difference over 
(dB) 0.1 2 


a 


ae 
209.94 | _¥es 9 
a _ a 
| vee 
| Xess 
| wee 
[1 [ure _ Youn 

__ ses 


Magnitude iors Sante 2dB are not listed. 





difference in dB among all the layers is listed if it is greater than 2. If the phases of 
the coefficients deviate more than 0.17 in any layer, that particular mode is also 
singled out. The location of the mode of the revised program is within a relative 
accuracy of 2~” achieved through first order Newton-Raphson iteration. Even 
though discrepancies still exist when the duct is 28 meters or higher, it is clear that 
the revised program computes more accurately than the original one. 

For the cases where the two methods of evaluating the A; and B; coefficients 
disagree, it has been observed that one of the methods always leads to A; values 
which are little changed when the accuracy in mode location is varied, while the 


other method produces A; values which shift toward the results of the other method 
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as the accuracy of mode location improves. Based on this observation, a consistency 
check is implemented into the program to identify the method which converges 
better. For the 120 cases investigated, when this method of faster convergence is 
used, the A, coefficients obtained from mode locations estimated with or without the 
Newton-Raphson first order iteration differ only by 0.06 dB in magnitude and 
0.00137 in phase at most. In fact, few cases show differences more than 0.002 dB 
and 0.00017. This allowed the Newton-Raphson iteration to be removed in this 
revision. 

Table 2 compares the performance between the original and the revised 
programs. The time spent to find the modes has been reduced by an average of 
22.58%. The revised program can always produce the modes found by the original 
program. Moreover, the mode search is stable for the new program: the time it 
requires to search for the modes is about the same for similar profiles. The sudden 
jumps in mode search time for the 24 m and the 40 m cases, which indicate troubles 
during the search, no longer happen. 

With the proper method of evaluating the A ; and B; coefficients determined by 
the consistency check, the output of the revised program differs from the original 
program in some cases. The most serious deviation has been observed for the 38 m 
duct height case as shown in Tables 3 and 4. For example, at a range of 36.5 km 
with the transmitter at a height of 25 m and the receiver at 10 m, the coherent path 
loss is 175.93 dB from the original program, and is 167.90 dB from the revised 


program. 
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TABLE 2 OVERALL MODE SEARCH PERFORMANCE COMPARISON 


DUCT ORIGINAL 
HEIGHT 


Time 
rovement 





fe 


© 
S 
oe) 
XQ 
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ae) 
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H 
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_ 
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12 1:09:40 86 1:01:44 | go | 11.39% 


6 

“s 
3:42:59 

28 2:07:42 
11 


a | 


Total 40:23:54 Biligt:22 22.58% 
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TABLE 3. ORIGINAL PROGRAM 38 M DUCT OUTPUT 


frequency = 9600.0000 mhz 
range zt ze coherent incoherent coherent incoherent horizon 
(km) (m) (m) mode sum mode sum path loss path loss (km) 
(dB ) 


rape 
25. 
Coe 
2 
aay 
Fa i 
23% 
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TABLE 4. REVISED PROGRAM 38 M DUCT OUTPUT 


frequency = 9600.0000 mhz 
range zt zr coherent incoherent coherent incoherent horizon 
(km) (m) (m) mode sum mode sum path loss path loss (km) 
(dB) (dB) (dB) (dB) 

25. =o 155. 156. 
me 140. 143. 
-4, 142. 145. 
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145. 
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B. Recommendations 

The mode search protocol of this program needs to be revised. Since the 
search is limited by the number of modes to be found and the maximum range 
attenuation rate accepted, it is more logical to begin with locating the mode of the 
lowest attenuation, and then proceed to look for the next one in the order of 
increasing attenuation rate. Furthermore, there appears to be only a single ‘phase 
line’ of vanishing real part of the modal function on which all the modes are located. 
This line extends from lower to higher range attenuation rates. The partition of the 
search region into rectangles, as has been done in this program, tends to cut the 
"phase line” into segments before the program starts to search for the end points of 
these segments and then follow the segments in different directions. It is clear that 
a better way for mode iach is to find the lower end of the single "phase line" then 


follow it to the other end. 
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APPENDIX A: SUBROUTINE XCADD 


This Appendix lists the addition subroutine XCADD which returns the complex 
exponent of the sum when the complex exponents of the addends are given. This is 


a complete re-write of the original subroutine of the same name. 
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subroutine xcadd(zx,z1x,Z2x) 


c 
fe Given 21x and zéx, this subroutine adds the two complex numbers 
Cc z1=exp(z1x) and z2=exp(z2x) for z=exp(zx) and returns 2x. 
c 
Cc inputs... 
Cc z1x=complex exponent of the complex number 21 
Cc z2x=complex exponent of the complex number 22 
Cc 
Cc outputs... 
ec zx=complex exponent of the complex number z 
c 
c subroutines called... 
c 
CRRRKKKKKEK 
implicit real*8 (a-h,o-z) 
complex*16 zx,z1x,z2x,zt1x,zt2x,clogzh,dsum,czero,cerrx,cone,chpi 
parameter (pi=3. 141592653589 793238462643d0, twopi=2.d0*pi, 
+ hpi=0.5d0*pi,zero=0.d0,c16=1.d0/6.d0, 
+ bit14=1.d0/16384.d0,bit24=bit14/1024.d0,ctol=bit14, 
+ = dpi=2259.d0/4294967296.d0/4294967296.d0, hdpi=dpi/2.d0, 
+ e2m54=-3.742994775023704819d1 , e2p27=-0.5d0*e2m54, 
+ chpi=(0.d0,1.57079632679489661923132d0), cone=(1.d0,0.d0), 
+  czero=(0.d0,0.d0),cerrx=(-3.742994775023704819d1,0.d0)) 
Cc cerrx=e2m54=-54*log(2)=exponent below machine accuracy 
dimension ztmpo(2),stmp(2) 
equivalence (ztmp,clogzh),(stmp,dsum) 
Ch Rae 
Cc Replace the input variables with a local variable so that 
Cc equations in the form of y=xt+y will not lead to confusion. 
Cc 
zt1x=z1x 
Zt 2x=Z2x 
Cc 
clogzh=0.5d0*(zt1x-zt2x) 
dxh=ztmp(1) 
if(dxh .lt. zero) then 
ZX=Zt2Xx 
dxh=-dxh 
else 
zx=zt1x 
end if 
c* kkkk&kkk 
Cc machine accuracy = 2**(-53) 


C 2** (27 )=e**e2p27 
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if (dxh .ge. e2p27) then 
return 
else 
2x=0.5d0*(zt1x+zt2x) 
dsunecdexp(clogzh) 
dsum=1.d0/dsum+dsum 
if (cdabs(dsum) .gt. ctol) then 
zx=cdlog(dsum)+zx 
else 
Cancellation is serious. Im{clogzh] is close to pi/2 or -pi/2. 
yi=dnint(ztmp(2)/twopi )*2.d0 
ztmp(2)=ztme(2)-pi*yi 
dyi=dpi*yi 
if (ztmp(2) .lt. zero) then 
clogzh=-clogzh 
dyi=-dyi 
end if 
ztmp(2)=(ztmp(2)-hpi )-hdpi-dyi 
dsum=2.d0*clogzh* (cone+c16*clogzh*clogzh) 
if (dsum .eq. czero) then 
Note that a complete cancellation of two nonzero numbers of 
order one is considered to be as accurate as what is allowed 
by the machine and the algorithm. 
zx=cerrx+chpi+zx 
else 
dsum=cdlog(dsum) 
if (stmp(1) .lt. e2m54) stmp(1)=e2m54 
zx=dsum chpi+zx 
end if 
end if 
return 
end if 


end 


4] 


APPENDIX B: SUBROUTINE FZEROX 


This Appendix includes the listing of the subroutine FZEROX which searches 
identifies the meshes which may contain modes within a contour rectangle. The 


Shellman-Morffit mode locating algorithm has been completely replaced. 
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Cc 
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subroutine fzerox(tleft,tright, tbot, ttop, tmsh0, zeros,ni,nf) 


Creeee 


fzerox is a routine for finding the zeroes of a complex function, f, 
which lie within a specified rectangular region of the 
complex qil plane, assuming that the function has only 


simple zeroes over this rectangle. 


parameters specifying the search rectangle: 

tleft - value of the real part of qiil at the left edge. 

tright- value of the real part of qi1 at the right edge. 

tbot - value of the imaginary part of qiil at the bottom edge. 
(this is set to 0.) 

ttop - value of the imaginary part of qil at the top edge. 

tmesh - set equal to about half the average spacing between 
zeroes within the rectangle. A smaller value may be used 
as a safety measure, but too small a value will result 
in excessively long run time. 

zeros - output list of (complex) values of qil at which 
zeroes are found. 


nf-ni - the number of zeroes found 


subroutines cal ledd-- 
f indfx 
roots 


nomshx 
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implicit double precision (a-h,o-z) 

complex*16 £10, f01,f11, fxnew, fxold, fx00, fx10, fx01, fx11, 
+ czero,one,ci,sol, zeros 

parameter (czero=(0.d0,0.d0), one=(1.d0,0.d0),c1=(0.d0,1.d0)) 


31 $include: ‘'mlaparm.inc' 


***** Begin listing of: mlaparm.inc 


include file to define the 
maximum # of layers (mxlayr) 
maximum # of modes (mxmode) 


parameter (mxlayr=35 ) 
parameter (mxmode=127) 


**¥*** End listing of: mlaparm. inc 


Me ce 
2c 
eee 
4 ¢ 
a. EC 
6 

7 
32 
33) ¢ 
34 c¢ 
255] 
36°C 


dimension kedge1(100),kedge2( 100), kedge3( 100), kedge4( 100), 

+  loctér¢(mxmode), loc12i (mxmode), loc23r (mxmode) , loc23i (mxmode), 
+ loc34r(mxmode), loc34i(mxmode), loc41r(mxmode), loc41i (mxmode), 
+  sol(3),theta(2),zeros(2*mxmodet+1) 
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Cxxxee 


chelatehehel 


Cxrene 


Cc 


Cc 


common /tmccom/tmesh 


maxnsq - maximum number of mesh squares allowed on any one 


phase line 


maxnt - maximum number of times fzerox will reduce tmesh 


maxnsq=3*max0(int((ttop-tbot)/tmsh0), int((tright-tleft)/tmsh0) ) 


maxnt=2 


tmesh = tmsh0O 


i) 
(=) 


ntime 


go to 7 
tmesh=tmesh/2.0d0 
ntime = ntime+1 


if(ntime .gt. maxnt) go to 97 


continue 


calculate coordinates of rectangle edges in tmesh units 


jlt = idnint(tleft/tmesh-0.5d0) 

jrt = idnint(tright/tmesh+0.5d0) 
Jtop = idnint(ttop/tmesh+1.5d0) 

jbot = 0 


initialize parameters for starting search at upper left 


corner of search rectangle 


ki = jtop 
kr = jlt 
kedge = 1 
call findfx(kr,ki, fxnew, xnew, ynew) 
nrei=0 
nre2=0 
nre3=0 
nre4=0 
knot 12=0 
knot23=0 
knot 34=0 
knot41=0 
nf=ni 


nit=ni+1 
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Cxexen 


10 


15 


20 


Cuxene 


21 


23 


Curae 


Cc 


go to 15 


continue 

1fcnrzt. te) go to 15 
write(16,2000) nrzlt 

go to 5 

nrzl=0 

nrsqu = 0 

f xold=fxnew 

xold=xnew 

yold=ynew 

go to (21,26,31,36),kedge 


search along left edge of rectangle for changes in the 


sign of imag(f) 


continue 
if(ki.eq.jbot) then 
kedge=2 
go to 26 
end if 
ki = ki-1 
call findfx(kr,ki, fxnew, xnew, ynew) 
if (yold*ynew .gt. 0.d0) go to 20 
if(nrel.eq.0) go to 23 


check if crossing point has been previously found 


do 22 k=1,nre1 
if(ki.eq.kedgel(k)) go to 20 


cont inue 


follow phase Line through rectangular region 


fx01=fxold 
fx01r=xold 
fx01i=yold 
f x00=f xnew 
fx00r=xnew 
f x001=ynew 
li = ki 

lr = jlt 
go to 43 


search along bottom edge of rectangle for changes in the 
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levee sign of imag(f) 


128 ¢ 

129 26 continue 

130 if¢kr.eq.jrt) then 

131 kedge=3 

132 go to 31 

133 end if 

134 kr = kr+1 

135 call findfx(kr,ki, fxnew, xnew, ynew) 
136 if (yold*ynew .gt. 0.d0) go to 20 
137 if(nre2.eq.0) go to 28 

138 c 

139 c¢ check if crossing point has been previously found 
140 c¢ 

141 do 27 k=1,nre2 

142 if (kr.eq.kedge2(k)) go to 20 

143 27 cont inue 

144 Cc 

145 c¢ follow phase line through rectangular region 
146 c¢ 

147 28 fx00=fxold 

148 fx00r=xold 

149 fx00i=yold 

150 fx10=fxnew 

151 fx10r=xnew 

152 fx10i=ynew 

153 li = jbot 

154 lr = kr-1 

155 go to 48 


156 ChNKKE 

157 c search along right edge of rectangle for sign changes in imag(f). 
1565°c 

159 31 continue 


160 if(ki.eq.jtop) then 

161 kedge=4 

162 go to 36 

163 end if 

164 ki = ki+1 

165 call findfx(kr,ki, fxnew, xnew, ynew) 
166 if (yold*ynew .gt. 0.d0) go to 20 
167 if(nre3.eq.0) go to 33 

168 Cc 

169 c¢ check if crossing point has been previously found 
N7O¢ 

171 do 32 k=1,nre3 
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172 if(ki.eq.kedge3(k)) go to 20 
173 32 continue 


174 Cc 

iy os “e follow phase line through rectangular region 
176 c¢ 

iv 335 fx10=fxold 
178 fx10r=xold 
179 fx10i=yold 
180 fx11=fxnew 
181 fx11r=xnew 
182 fx11i=ynew 
183 li = ki-1 
184 lr = jrt-1 
185 go to 53 


186 cxxxxx 

187 c search along top edge of rectangle for sign changes in imag(f). 
188 c¢ 

189 36 continue 


190 if(kr.eq.jlt) go to 80 

191 kr = kr-4 

192 call findfx(kr,ki, fxnew, xnew, ynew) 
193 if (yold*ynew .gt. 0.d0) go to 20 
194 if(nre4.eq.0) go to 38 

vo: Cc 

196 Cc check if crossing point has been previously found 
ise c 

198 do 37 k=1,nre4 

199 if(kr.eq.kedge4(k)) go to 20 

200 37 cont inue 

201 ¢ 

202 c follow phase line through rectangular region 
203 c¢ 

204 38 fx11=fxold 

205 fx1lir=xold 

206 fx1liz=yold 

207 fx01=fxnew 

208 fx01r=xnew 

209 fx01i=ynew 

210 li = jtop-1 

211 lr = kr 

212 go to 58 


213 Cree 

214 c enter mesh square from left side or exit rectangle at right edge. 
215 

216 41 tr=(r+1 
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217 1f Clr vlew jrt-1) go to 42 


218 nre3=nre3+1 

219 kedge3(nre3)=li+1 

220 go to 10 

221 42 fx01=fx11 

222 fxO1r=fxiir 

2235 fx01i=fx111 

224 f x00=fx10 

225 fx00r=fx10r 

226 fx00i=fx101 

227 43 continue 

228 cell findfx(lr+1,li+1, fx11, fxiir, fx111) 

229 call findfx¢(lr+1, li, fx10, fx10r, fx101) 

230 CrKKKKKK 

231 Cc Determine the edge of exit of im(f)=0 from current mesh. 
232 edgeit=fx01i*fx111 

233 edgeib=fx00i*fx10i 

234 if Cedgeib .gt. 0.d0) then 

255 0c Im(f)}=0 goes through the 01 to 10 line. 

236 if (edgeit .gt. 0.d0) then 

25a °c Im(f)}=0 goes through the 10 to 11 edge (edge 1). 
238 lout=1 

239 else 

240 9c Im(f)=0 goes through the 01 to 11 edge (edge 2) 
241 lout=2 

242 end if 

243 else 

244 Cc Im(f}=0 goes through the 00 to 10 edge (edge 4) 
245 lout=4 

246 if (edgeit .lt. 0.d0) then 

C47 6c Im(f)}=0 also runs through 01 to 11 and 10 to 11 edges. 
248 ¢ Store crossing location and in/out information. 
249 knot 34=knot34+1 

Zo0.c loc34r(knot34)=lr 

231 ne loc341(knot34)=Li 

252 end if 

2o5 end if 

254 kk kik 

255 go to 60 


256 Cxrene 
257 c enter mesh square from bottom side or exit rectangle at top edge. 
258 46 Li=li+1 


259 if (li .le. jtop-1) go to 47 
260 nre4=nre4+1 
261 kedge4(nre4)=lr 
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262 
263 
264 
265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
co 
276 
277 
278 
279 
280 
281 
282 
283 
284 
285 
286 
287 
288 
289 
290 
291 
292 
293 
294 
295 
296 
297 
298 
299 
300 
301 
302 
303 
304 
305 
306 


go to 10 
47 f x00=f x01 
fx00r=fxOir 
fx00i=fx011 
fx10=fx1i1 
fx10r=fxiir 
fxi0i=fxiti 
48 continue 
call findfxclr,lit+i, fx01, fxOir, fx011) 
call findfx(tr+1, litt, fx11, fxiir, fxili) 
CXREKKKK 
c Determine the edge of exit of im(f)=0 from current mesh. 
edgei l=fx00i*fx01i 
edgeir=fxi0i*fx11i 
if Cedgeir .gt. 0.d0) then 


Cc ImCf)=0 goes through the 00 to 11 line. 
if (edgeil .gt. 0.d0) then 
Cc Im(f)=0 goes through the 01 to 11 edge (edge 2) 
Lout=2 
else 
Cc Im(f)=0 goes through the 00 to 01 edge (edge 3). 
Lout=3 
end if 
else 
Cc Im(f)=0 goes through the 10 to 11 edge (edge 1) 
lout=1 
if Cedgeil .lt. 0.d0) then 
c Im(f)=0 also runs through 00 to 01 and 01 to 11 edges. 
c Store crossing location and in/out information. 
knot41=knot41+1 
Cc Loc4ir(knot41)=lr 
(© Loc41i(knot4i)=Li 
end if 
end if 
CXRKKKRK 
go to 60 
Ck ke 


c enter mesh square from right side or exit rectangle at left edge. 


51 lr=lr-1 
1f Clr «ge. slit) go to 52 
nrei=nrei+1 
kedgel(nrei)=li 
go to 10 

52 fxi1=fx01 
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307 
308 
309 
310 
311 
312 
313 
314 
315 
316 
317 
318 
319 
320 
321 
322 
323 
324 
325 
326 
327 
328 
329 
330 
331 
332 
333 
334 
335 
336 
337 
338 
339 
340 
341 
342 
343 
344 
345 
346 
347 
348 
349 
350 
351 


fx11ir=fx0ir 
fxi1i=fx01i 
fx10=fx00 
fx10r=fx00r 
fx10i=fx00i 
53 continue 
call findfx(lr,li+1, x01, fxO1r, fx011) 
call findfx(lr,l1,fx00, fx00r, fx001) 


Cxxeaknn 
Cc Determine the edge of exit of im¢f)=0 from current mesh. 
edgeit=fx01i*fx11i 
edgeib=fx001*fx101 
if (edgeit .gt. 0.d0) then 
Cc Im(f)=0 goes through the 01 to 10 line. 
if (edgeib .gt. 0.d0) then 
c Im(f)=0 goes through the 00 to 01 edge (edge 3). 
lout=3 
else 
c Im(f)=0 goes through the 00 to 10 edge (edge 4) 
lout=4 
end if 
else 
c Im(f)=0 goes through the 01 to 11 edge (edge 2) 
lout=2 
if (edgeib .lt. 0.d0) then 
Cc Im(f)=0 also runs through 00 to 10 and 00 to 01 edges. 
Cc Store crossing location and in/out information. 
knot 12=knot 12+ 1 
C locilér¢knoti2)=lr 
c loc12i(knot12)=Li 
end if 
end if 
Ck iO eek 
go to 60 
Co ike 


c enter mesh square from top side or exit rectangle at bottom edge. 
56 Li=li-1 

if (li .ge. jbot) go to 57 

nre2=nre2+1 

kedge2(nre2)=Ir+1 

go to 10 
57 fx01=fx00 

fx01r=fx00r 

fx01i1=fx001 

fx11=fx10 
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352 
353 
354 
355 
356 
357 
358 
359 
360 
361 
362 
363 
364 
365 
366 
367 
368 
369 
370 
371 
372 
373 
374 
375 
376 
377 
378 
379 
380 
381 
382 
383 
384 
385 
386 
387 
388 
389 
390 
391 
392 
393 
394 
395 
396 


fx11r=fx10r 
fx11i1=fx101 
58 continue 
call findfx(lr, li, fx00, fx00r, fx001) 
call findfx(lr+1, li, fx10,fx10r, fx101) 
CR kee 
Cc Determine the edge of exit of im(f)=0 from current mesh. 
edgeil=fx00i1*fx011 
edgeir=fx10i*fx111 
if (Cedgeil .gt. 0.d0) then 
Cc Im(f)=0 goes through the 00 to 11 line. 
if (edgeir .gt. 0.d0) then 
Cc Im(f)=0 goes through the 00 to 10 edge (edge 4) 
lout=4 
else 
Cc Im(f)=0 goes through the 10 to 11 edge (edge 1). 
Lout=1 
end if 
else 
Cc Im(f)=0 goes through the 00 to 01 edge (edge 3) 
lout=3 
if Cedgeir .lt. 0.d0) then 
Cc Im(f)=0 also runs through 00 to 10 and 10 to 11 edges. 
Cc Store crossing location and in/out information. 
knot23=knot23+1 
c loc23r(knot23)=lr 
c loc23i(knot23)=1li 
end if 
end if 
c 
CR IR Ree 
60 continue 
nrsqu=nrsqu+ 1 
if(nmrsqu .gt. maxnsq) go to 95 
Co Ra ie 
c Test for there being at least one re(f)=0 line entering and 
Cc leaving the mesh square. 
c 
17 (CCPxODe*{x10r ..9gt-.0.d0) Vand. sCix0ir*{xiir gt. 0.00) 
+ and. (fx00r*fx0ir .gt. 0.d0)) go to (41,46,51,56) lout 
c 
c Computate the values of the modal function at the corners of a 
Cc a mesh square to determine its Taylor series to the 3rd order 
Cc for estimating its root locations. 
c 


397 
398 
399 
400 
401 
402 
403 
404 
405 
406 
407 
408 
409 
410 
411 
412 
413 
414 
415 
416 
417 
418 
419 
420 
421 
422 
423 
424 
425 
426 
427 
428 
429 
430 
431 
432 
433 
434 
435 
436 
437 
438 
439 
440 
441 


c f00=one 

f 10=cdexp( fx10- fx00)-one 

f01=cdexp( fx01-fx00)-one 

f 11=cdexp( fx11-fx00)-one 
c 
CER KEKRKRE REE EKEEKEE KEE EE KERR EERE EREKEKEREEEEEEEEEEEREREEKEE 
Cc write (16,3001) ni,nf,t(r,li,knot12,knot23, knot34,knot41 
c 3001 format(/' ni, nf, tr, li and knot12, 23, 34 and 43 before ROOTS 
c + 3'/, 216,2x,216,2x,416) 
c 
crxseeeeeeex estimate locations of zeroes by radicals “"=*777tr st -=7578 
c 

call roots(f10,f01,f11,sol,nrsol ) 


do 63 n=1,nrsol 
ureal = dreal(sol(n)) 


dimag(sol(n)) 


ulmag 
if cureal .lt. 0.d0 .or. ureal .gt. 1.0d0) go to 63 
if (uimag .lt. 0.d0 .or. uimag .gt. 1.0d0) go to 63 
62 theta(1)=(lr+ureal )*tmesh 
theta(2)=(lit+uimag)*tmesh 
nf = nf+1 
zeros(nf )=demp|x(theta(1), theta(2)) 
nrzl=nrzl+1 - 
63 cont inue 
CREEK KEKKEKEKK KK KEKE KEEKKEKEKEKEKKEKEKKEEKEKKEKKEEEKEKKKKKKE 
c write (16,3002) ni,nf,nrsol 
c 3002 format(/' out of ROOTS at 63, ni, nf and # of roots ',3i14) 


CREAR KEKKEKKKKKE KEKE KKK KKK EKER KEEKEKEKKEEKKKKKKKEKKKEKEKKEKKKEK 


c continue following the phase line 
go to (41,46,51,56) lout 

(chghelelehelel 

cc 

80 cont inue 

c 
return 

Ccxxnnn 

95 cont inue 
write(16,9500) 


write(16,4001)Ir,li,ni,nf, tmesh 
write(* ,9500) 
4001 format('go to 5 from 95 at tr, ti =",16,',"716)" mi, nr = 16, 
+'.',16,', mesh size =',d14.6) 


go to 5 


Curnne 


DZ 


442 97 continue 


443 write(16,9700) 

444 write(16,4002)\r,li,ni,nf,tmesh 

445 write(* ,9700) 

446 4002 format('go to 5 from 97 at lr, li =',16,',',i16,' ni, nf =',i6, 
447 +','16,', mesh size =',d14.6,/'zeroes found are kept.') 
448 ¢ nf=ni 

449 c¢ 

450 return 

451 Cc 

452 c**** format statements 

453 9500 format(/5x,'too many squares on same phase line -- ', 
454 $ ‘reduce tmesh and start over') 


455 9700 format(/5x,'tmesh has been reduced but problems remain in', 
456 $ ' executing fzerox') 

457 c¢ 

458 end 


By. 


APPENDIX C: SUBROUTINE ROOTS 


This Appendix contains the listing of the subroutine ROOTS. This subroutine 
replaces the portion of the subroutine FZEROX where the coefficients of a quadratic 
equation are determined, and the subroutine QUAD for locating the zeroes of a 
quadratic polynomial. In the revised subroutine FZEROX, the roots of a cubic 


polynomia! has to be found. This subroutine determines these zeroes by radicals. 
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subroutine roots (f1,f2,f3,sol,nrsol) 
CREAR KKAEKREKREKEKEEKEREK KE EKEEEKEEEEKEKEEKKEEEKEKEKREKRKEKEEKEEEREEEREKE 
This subroutine finds the roots of a third order polynomial by 
radicals when the values of this polynomial at z=0, z=1, z=i and 
z=1+i are given as f0=1, f1+f0, f2+f0 and f3+f0 respectively. 
Note that this algorithm takes cubic roots of two complex numbers 
(hence the name ‘solution by radicals') and use their linear 


a a a aaa 


combinations as the roots of a third order polynomial. 


CREAR RKE KEKE KEKE KKK AKER EKER EKER KEREEREREEREKREEKKKKKK 
implicit real*8 (a-h, o-z) 
complex*16 f1,f2,f3,zero,one,ci,epi14,em14, ep23, em23, 


+ fa, fb,fc,fd,fal, fae, fa3, fais,p,q,delt,z,zm,u,v,sol 
parameter (xbit52=52.d0*0.69314718055994531d0, thrd=1.d0/3.d0, 

+ bit50=1.d0/33554432.d0/33554432.d0,bit51=bit50/2.d0, 

oa bit52=bit51/2.d0, tol=0.001d0, 

+ zero=(0.d0,0.d0), one=(1.d0,0.d0),ci=(0.d0,1.d0), 

+ epi4=(0.5d0,0.5d0) ,em14=(0.5d0,-0.5d0), 

+ ep23=(-0.5d0,0.86602540378443864675d0), 

+ em23=(-0.5d0, -0.86602540378443864675d0) ) 


dimension sol(*) 
fa=one 
fb=( f2-ci*fit+emi4*f3) 
fc=((epi4tone)*f1- (em14+one)* f2+ci* f3) 
fd=(emi4*( f2-f1)-ep14* f3) 
if (cdabs(fb) .le. bit50) fb=zero 
if (cdabs(fc) .le. bit51) fc=zero 
if (cdabs(fd) .le. bit52) fd=zero 
if (fd .ne. zero) then 
fai=(-thrd)*fc/fd 
fa2=fb/fd 
fa3=fa/fd 
fais=fai*fal 
p=thrd* fa2- fais 
q=0.5d0*( fa3+fal*fa2)-fal*fals 
if (p .eq. zero) then 
if (q. eq. zero) then 
nrsol=1 
sol(1)=fal 
return 
else 
nrsol=3 
u=((-2.d0)*q)**thrd 
sol(1)=u+fal 
sol(2)=ep23*ut fat 
sol(3)=em23*utfal 
return 
end if 
else 


=) 


if (q. eq. zero) then 
nrsol=3 
sol(1)=fal 
u=cdsqrt((-3.d0)*p) 
sol(2)=faitu 
sol(3)=fai-u 
return 
else 
v=p/q 
z=p*v*v 
abs z=cdabs(z) 
if (absz .lt. tol) then 
zm=-Z 
fn=dint(1.d0-xbit52/dlog(absz) ) 
lastn=idint(fn)-1 
dnn=fn-0.5d0 
dnd=fn+1.0d0 
delt=one 
do 100 nt=1,lastn 
dnn=dnn- 1.d0 
dnd=dnd- 1.d0 
del t=(dnn/dnd)*del t* zmtone 
100 cont inue 
del t=(0.5d0*delt/q)**thrd 
u=p*delt 
v=-1.d0/delt 
else 
del t=cdsqrt (one+z)-one 
u=(q*delt)**thrd 
v=-p/u 
end if 
nrsol=3 
sol(1)=u+v+ fal 
sol(2)=ep23*utem23*v+ fal 
sol(3)=em23*ut+ep23* vtfal 
return 
end if 
end if 
else if (fc .me. zero) then 
if (fb .eq. zero) then 
if (fa .eq. zero) then 
nrsol=1 
sol(1)=zero 
return 
else 
nrsol=2 
z=cdsaqrt(-fa/fc) 
sol(1)=z 
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200 


sol(2)=-z 
return 
end if 
else 
fa1=0.5d0* fb/ fc 
fa2=fa/fc 
z=fa2/fai/fal 
abs z=cdabs(z) 
if (absz .lt. tol) then 


fn=dint(1.d0-xbit52/dlog(absz)) 


lastn=idint(fn)-1 

dnn=fn-0.5d0 

dnd=fn+1.0d0 

del t=one 

do 200 nt=1,lastn 
dnn=dnn- 1.d0 
dnd=dnd- 1.d0 


del t=(dnn/dnd)*delt*z+one 


continue 
delt=-0.5d0*delt/fal 
nrsol=2 
sol(1)=fa2*delt 
sol(2)=1.d0/delt 
return 

else 
del t=cdsqrt(one-z) 
nrsol=2 


sol(1)=-fal*(one-delt) 
sol(2)=- fal* (one+delt) 


return 
end if 

end if 

else if (fb .ne. zero) then 
nrsol=1 
sol(1)=-fa/fb 
return 

else 
nrsol=1 
sol(1)=ep14 
return 

end if 

end 


ay 


APPENDIX D: SUBROUTINE ABCOEF 


This Appendix contains the listing of the subroutine ABCOEF. The consistency 
self-checking procedure has been implemented to determine the correct method to 


evaluate the A; and B; coefficients. 
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subroutine abcoef(zero,m) 


CR Rk 
c For each mode m, this suboutine calculates A-B coefficients in 
C all layers for combining two linearly independent solutions of 
c Stokes' equation to form the height gain function: 

C 

Cc height gain=exp(becoefx(l,m))*(k1*exp(acoefx (l,m) )+k2) 

c 

Cc where k1 and k2 are two independent solutions to Stokes’ 

C equation. In the top layer (i.e. nzlayr) the height gain is: 
G 

Cc height gain=exp(beoefx(l,m))*h2 

le where h2 is a solution to the Stokes' equation associated 

Cc with outgoing energy flow. Here ki and k2 are proportional 

c to the k1 and k2 used by Marcus and the h2 is proportional 

Cc to a modified Hankle function of order 1/3. 

G Inputs... 

Cc zero-an eigenvalue in qi1 space 

Cc outputs... 

G acoefx-two dimensional array of complex exponents 

c coefficients used to combine two linearly 

Cc independent solutions of stokes' equation 

(e bcoefx-two dimensional array of complex exponents 

c coefficients used for normalizing the height gains 
C note: acoefx and bcoefx are passed by the 

c common block /pap2/ 

Cc subroutines called... 

c xcdal 

Cc xcadd 

Cc common block areas... 

C com) 

Cc com2 

c pap1 

fe pape 

Cee Ree 


implicit real*8(a-h,o-z) 
complex*16 acoefx,bcoefx,cqij,h2xq1,dh2xq1,h2xq2, dh2xq2,k1xq1, 
$ dkixql,k1ixq2,dk1xq2,k2xq1,dk2xq1,k2xq2, dk2xq2,h2dk1x, 
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46 $ dh2k1x,h2dk2x,dh2k2x , numax, denax, numbx, denbx, int 1x, int2x, 

47 $ hyx,dhyx,kidhyx, dk ihyx,dk2hyx,k2dhyx, gamma, dgamdq,i, 

48 $ koal23,rtsum,zero,ql,q2,sum,surfno,dqij,dqijdz,sqng, 

49 $ dnumbx,dhux,dhlx,e13x,cneg, cldqzl,cldqzm,cigama, koawav, tthd, 
50 + tacoef,dacoef 

51 

52 parameter (downi=1.d-3,downr=1.d-3/0.43429448190325 18d0, 

53 + pi=3.141592653589793238462643d0, 

54 + 1=(0.0d0,1.0d0),tthd=(2.d0/3.d0)*i, 

55 + cneg=(0.0d0,3.141592653589793238462643d0) , e13x=cneg/3.d0) 
56 (chehtehehakel 

SY iauae ° mxlayr=maximum number of layers allowed 

58 ¢ mxmode=maximum number of modes allowed 

59 

60 c¢ 

61 Cc use include file for parameters of 

62 Cc use include file for parameters of 

63°56 mxlayr max # layers 

64 ¢ mxmode max # modes 

657e 


66 include: 'mlaparm.inc'! 


*xe** Begin Listing of: mlaparm.inc 


ic 

Ze include file’to define the 

Lae © maximum # of layers (mxlayr) 

4 c¢ maximum # of modes (mxmode) 

aC 

6 parameter (mxlayr=35 ) 

if parameter (mxmode=127) 

*eee* End Listing of: mlaparm. inc 
67 ¢ 
68 c¢ 
69 ckeRRS 
(Ady fe acoefx-two dimensional complex array used for combining two 
71 Cc independent solutions to stokes' equation 
l2 “ec bcoefx-two dimensional complex array used for normalizing height 
73 Cc gain 
74 =O Cqij-two dimensional array containing coefficients for evaluating 
15 € qij in terms of qil 
76 C¢ dqij-array containing coefficients for evaluating qij in terms of 
77 oc qii 
78 dqijdz-array containing derivatives of qi(z) in the different 
79 Cc layers 
80°C Zi-array containing input hesights for the modified refractivity 
81 
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dimension acoefx(mxlayr,mxmode), 
$ bcoefx(mxlayr ,mxmode), 
$ dqij(mxlayr),cqij(mxlayr,2),dqijdz(mxlayr),zi(mxlayr+1) 


heltehelehel 


common /comi/freq,waveno,sqng 

common /com2/cqij,dqij,dqijdz,nzlayr 
common /pap1/nrmode,koa123,surfno,zi 
common /pap2/acoefx, bcoef x 


Cuxeee 


Cc 


check for single layer 


set a complex variable koawav=-i*koal23/(waveno*waveno) to 


avoid repeating computations 


koawav=- 1*koa123/(waveno* waveno) 


if(nzlayr .eq. 1)then 
qil=cqij(1,1)+zero*dqij(1) 
call surf(qi,gamma,dgamdq) 
call xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 
dh2xq1=dh2xq1i+e13x 
int 1x=cdlog(koawav*dgamdq-q1/dqi jdz(1))+2.0d0*h2xq1 
int2x=2.0d0*dh2xq1-cdlog(-dai jdz(1)) 
call xcadd(sumx, int1x, int2x) 
rtsumx=0.5d0*sumx 
beoefx(1,m)=-rtsumx 
return 
end if 


cldqzl=cdlog(-dqi jdz(1)) 


if | equals one then initialize cumulants and caculate a's and 


b's in bottom layer using ground boundary conditions. 


qi=cqij(1,1)+zero*dqi j(1) 

call xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 
dk2xq1=dk2xqi+cneg 

dk 1xq1=dk 1xq1-e13x 

call surf(q1, gamma, dgamdq) 

cigama=cdlog(i*gamma) 

call xcadd(numax, cldqzl-cneg+dk2xq1,cigama+cneg+k2xq1) 
call xcadd(denax,cigamat+k1xq1,cldqzl+dk1xq1) 


acoefx(1,m)=numax-denax 
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call xcadd(denbx,k2xq1,acoefx(1,m)+k1xq1) 
bcoefx(1,m)=-denbx 


calculate contributions to normalizing integrals. 


call xcadd(hyx,k2xq1, acoefx(1,m)+k1xq1) 
hyx=beoefx(1,m)+hyx 

call xcadd(dhyx,dk2xq1, acoefx(1,m)+dk1xq1) 

dh yx=be oef x(1,m)+dhyx 

int 1x=cedlog( koawav*dgamdq- q1/dqi jdz(1))+2.0d0*hyx 
int2x=2.0d0*dhyx-cldqzl 


call xcadd(sum, int1x, int2x) 


do 9010 t=2,nzlayr-1 


(m1i=l-1 

cldqzl=cdlog(-dqi jdz(l)) 

cldqzm=cdlog(dqi jdz¢(lm1)) 

qi=cqij(l,1)+zero*dgqij(l) 

call xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 
dk2xqi=dk2xqi+cneg 

dk 1xq1=dk1xq1-e13x 

q2=cqi j(im1,2)+zero*daqi j(im1) 

call xcdai(-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2 ,dh2xq2) 
dk2éxq2=dk2xq2+cneg 

dk 1xq2=dk1xq2-e13x 

call xcadd(hyx,k2xq2,acoefx(lm1,m)+k1xq2) 

call xcadd(dhyx,dk2xq2, acoefx(lm1,m)+dk1xq2) 

k idhyx=k1xq1+dhyx 

dk Ihyx=dk Ixqit+hyx 

dk2hyx=dk2xq1+hyx 

k2dh yx=k2xq1+dhyx 

call xcadd(denax, cldqzmkidhyx, cldqzl+dkthyx) 

call xcadd(numax, cldqzl-cneg+dk2hyx, cldqzm+cnegt+k2dhyx) 
acoef x(t ,m)=numax-denax 

call xcadd(denbx,k2xq1,acoef x(l ,m)+k1xq1) 
numbx=be oef x( lm1,m)+hyx 

dnumbx=becoefx(lm1,m)+dhyx 

beoefx(l ,m)=numbx-denbx 


calculate contribution to normalizing integrals. 


int1x=cdlog(-q1/dqi jdz(l)+q2/dqi jdz(lm1))+2.0d0*numbx 
call xcadd(sum,sumx, int1x) 

call xcadd(dhux ,dk2xq1,acoefx(l,m)+dk1xq1) 
dhux=be oef x( l ,m)+dhux 
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172 int 1x=2.0d0*dnumbx-cldqzm 


173 int2x=2.0d0*dhux-cldqzl 

174 call xcadd(sumx, sumx, int 1x) 

l75 call xcadd(sumx,sumx, int2x) 

176 9010 cont inue 

177 

178 

179 ¢ if | equals nzlayer, calculate a's and b's using outgoing 
180 c wave in top layer. 

181 

182 nzmi=nzlayr-1 

183 qi=cqi j(nzlayr, 1)+zero*dqi j(nzlayr) 

184 call xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 
185 dh2xq1l=dh2xq1+e13x 

186 q2=cqi j(nzmi, 2)+zero*dqi j(nzm1) 

187 call xcdai(-q2,k2xq2,dk2xq2,k1xq2,dk1xq2,h2xq2,dh2xq2) 
188 dk2xq2=dk2xq2+cneg 

189 dk 1xq2=dk1xq2-e13x 

190 call xcadd(hyx, kexq2, acoefx(nzm1,m)+k1xq2) 

191 numbx=bc oe fx(nzlayr-1,m)+hyx 

192 beoefx(nzlayr,m)=numbx-h2xq1 

193 

194 c¢ calculate contribution to cumulants. 

195 

196 int 1x=cdlog(-qli/dqi jdz(nzlayr )+q2/dqi jdz(nzm1) )+ 

197 $ 2 .0d0* numbx 

198 call xcadd(sum,sumx, int 1x) 

199 call xcadd(dhyx,dk2xq2,acoefx(nzm1,m)+dk 1xq2) 

200 dnumbx=bcoe fx (nzm1 ,m)+dhyx 

201 int 1x=2.0d0*dnumbx-cdlog(dqi jdz(nzm1)) 

202 call xcadd(sumx, sumx, int1x) 

203 dhux=be oe fx(nzlayr ,m)+dh2xq1 

204 Int2x=2.0d0*dhux-cdlog(-dqi jdz(nzlayr)) 

205 call xcadd(sum, sumx, 1nt2x) 

206 

207 c renormalize b's so that height gain integral equals unity. 
208 

209 rtsum=.5d0* sumx 

210 do 9000 (l=1,nzlayr 

211 beoefx(ll,m)=bcoefx(ll,m)-rtsumx 

212 9000 continue 

213 

214 CREAR ERE EERE EKER REE E REE KEKE EEE EEE 
215 

216 
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217 
218 
219 
220 
221 
222 
223 
224 
225 
226 
227 
228 
229 
230 
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 
248 
249 
250 
251 
252 
253 
254 
255 
256 
Zot 
258 
259 
260 
261 


l=nzlayr 

Imi=l-1 

cldqzm=cdlog(dqi jdz(lm1)) 
cldqzl=cdlog(-dqi jdz¢l)) 


c calculate q and associated quantities at bottom of layer l 


ql=cqij(l,1)+zero*dgi j(l) 
call xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1 ,dh2xq1) 
dh2xq1=dh2xqit+e13x 


q2=cqij(lm1,2)+zero*dgqi j(lm1) 

call xcdai(-q2,k2xq2,dk2xq2 ,k1xq2,dk 1xq2 ,hexq2, dh2xq2) 
dk2xq2=dk2xq2+cneg 

dk 1xq2=dk1xq2-e13x 


chahelehelel 


c Caculate acoefx(lm1,m),bcoefx¢(lm1,m) 
€ and cumulants using outgoing wave in nzlayr 
cerns 

dh2k 1x=dh2xqi+k1xq2 

h2dk1x=h2xq1+dk1xq2 

h2dk2x=h2xq1+dkexq2 

dh2k2x=dhexq1+kexq2 


call xcadd(denax,cldqzl-cneg+dh2k1x, cldqzm+cnegth2dk1x) 
call xcadd(numax, cldqzm+h2dk2x ,cldqzl+dh2k2x) 


c If in the nzlayr-1 layer the magnitudes of A coefficients from 
c integration up and down differ by less than 0.02 dB and their 
c phases differ by less than 0.001pi, the A and B coefficients 

c obtained from integration up will be accepted. 


tacoef=numax-denax 
dacoef=tacoef-acoefx(lm1,m) 
dif r=dabs(dreal (dacoef )) 

if (difr .lt. downr) then 
di f i=dimag(dacoef )/pi 
dif i=dabs(dif i-dnint(difi/2.d0)*2.d0) 
if (difi .lt. downi) return 

end if 


acoefx(lm1,m)=tacoef 


call xcadd(denbx ,k2xq2,acoefx(lm1,m)+k1xq2) 
beoefx(lm1,m)=h2xq1-denbx 
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Cc 


c 


cxrexe 


calculate contributions to cumulants 


sumx=cdlog(-qi/dqi jdz(l)+q2/dqi jdz(lm1))+2.0d0*h2xq1 
call xcadd(dhlx,dk2xq2, acoefx( lm1,m)+dk1xqQ2) 

dhl x=bcoe fx¢(lm1,m)+dhlx 

int 1x=2.0d0*dh2xq1-cldqzl 

call xcadd(int1x, sum, int1x) 

int2x=2.0d0*dhl x-cldqzm 

call xcadd(sum, int1x, int2x) 


do 9030 l=nzlayr-1,2,-1 
im1=l-1 

cldqzl=cdlog(-dqi jdzc(l)) 
cldqzm=cdlog(dqi jdz(lm1)) 


calculate q and associated quantities at bottom of layer lt 


ql=cqij(l,1)+zero*dqij cl) 

call xcdai(-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1,dh2xq1) 
dk2xql=dk2xq1+cneg 

dk 1xql=dk 1xq1-e13x 


q2=cqij(lm1,2)+zero*daqi j( lm1) 

call xcdai(-q2,k2xq2, dk2xq2,k1xq2,ak1xq2,h2xq2,dh2xq2) 
dk 2xq2=dk2xq2+cneg 

dk 1xq2=dk1xq2-e13x 

dh2xq2=dh2xq2+e13x 


Calculate acoefx(lm1,m),bcoefx(lm1,m) and cumulants 
using continuity relations in terms of the linearly 
independent functions k1 and k2 


call xcadd(hyx,k2xq1,acoefx(l,m)+k1xq1) 
call xcadd(dhyx,dk2xq1,acoefx(l,m)+dk1xq1) 
k Idhyx=k 1xq2+dhyx 

dk Ihyx=dk 1xq2+hyx 

Ak 2hyx=dk2xq2+hyx 

k2édhyx=k2xq2+dhyx 


call xcadd(denax,cldqzl-cneg+tk Idhyx, cldqzm+cneg+dk Thyx) 
call xcadd(numax,cldqzm+dk2hyx, cldqzl+k2dhyx) 
acoe fx( lm1,m)=numax- denax 


call xcadd(denbx,k2xq2,acoefx(lm1,m)+k1xqQ2) 


65 


307 numbx=be oe f x(t ,m)+hyx 


308 dnumbx=becoe f x¢ | ,m)+dhyx 

309 bcoe fx¢lm1 ,m)=numbx- denbx 

310 

311 ¢ calculate contributions to cumulants. 

312 

313 int 1x=cdlog(-qi/dai jdz(l)+q2/dqi jdz( lm1) )+2.0d0*numbx 
314 call xcadd(sumx,Sumx, int1x) 

315 call xcadd(dhlx, dk2xq2,acoefx(lm1,m)+dk1xq2) 

316 dhl x=beoefx(lm1,m)+dhlx 

317 int 1x=2 .0d0*dnumbx-cldqzl 

318 Int2x=2.0d0*dhl x-cldqzm 

319 call xcadd(sumx, sum, int1x) 

320 call xcadd(sumx,sumx, int2x) 

321 

322 9030 cont inue 

323 

324 cxreae 

B2> Cc if | equal to one calculate ground 

326 ¢ contribution to cumulants and renormalize bcoefx's 
327 

328 l=1 

329 ql=cqij(l,1)+zero*dgij(l) 

330 call xedai (-q1,k2xq1,dk2xq1,k1xq1,dk1xq1,h2xq1, dh2xq1) 
331 dkexql=dk2xqitcneg 

532 dk 1xq1=dk1xq1-e13x 

333 

334 call xcadd(hyx,k2xq1,acoefx(l,m)+k1xq1) 

335 call xcadd(dhyx,dk2xq1,acoefx(l,m)+dk1xq1) 

336 call surf(q1,gammna, dgamdq) 

337 numbx=becoe fx(l,m)+hyx 

338 dnumbx=bcoefx(l,m)+dhyx 

339 int1x=cdl og(koawav* dgamdq-q1/dqi jdz(l) )+2.0d0*numbx 
340 int2x=2.0d0*dnumbx- cdl og(-dqi jdz(1)) 

341 call xcadd(sumx, Sumx, 1nt1x) 

342 call xcadd(sumx,sumx, int2x) 

343 

344 ¢ renormalize b's so that height gain integrals equal unity. 
345 

346 rtsumx=.5d0*sumx 

347 

348 do 9020 Ll=1,nzlayr-1 

349 beoefx(ll,m)=bcoefx(ll,m)-rtsumx 

350 9020 cont inue 

351 
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352 
353 
354 
355 
356 


return 
end 


bcoe fx(nzlayr,m)=-rtsumx 
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